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Abstract In previous work, the numerical solution of the linearized gravitational field 
equations near space-like and null-infinity was discussed in the form of the spin-2 
zero-rest-mass equation for the perturbations of the conformal Weyl curvature. The 
motivation was to study the behavior of the field and properties of the numerical 
evolution of the system near infinity using Friedrich's conformal representation of 
space-like infinity as a cylinder. It has been pointed out by H.O. Kreiss and others that 
the numerical evolution of a system using second order wave equations has several 
advantages compared to a system of first order equations. Therefore, in the present 
paper we derive a system of second order wave equations and prove that the solution 
spaces of the two systems are the same if appropriate initial and boundary data are 
given. We study the properties of this system of coupled wave equations in the same 
geometric setting and discuss the differences between the two approaches. 



1 Introduction 

In a recent paper JTJ we studied the behaviour of gravitational perturbations on 
Minkowski space near space-like infinity. For this purpose we used the spin-2 zero- 
rest-mass equation which governs the perturbations of the Weyl curvature. The system 
was studied using a conformal gauge given by Friedrich |6| which represents space- 
like infinity as a cylinder connecting past and future null-infinity. This gauge exhibits 
explicitly the main issues of the propagation of gravitational waves in asymptotically 
flat space-times, namely the question of how is related to , or, in other words, 
how the late ingoing gravitational radiation affects the early outgoing radiation. 
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In JTJ it was shown how the field propagates near the cylinder. It turns out that 
the cylinder is a total characteristic in the sense that the equations turn into intrinsic 
equations on it, there appear no outward derivatives. Hence, the field on the cylinder 
is entirely determined by initial data on a space-like initial hyper-surface and there are 
no free data to be prescribed on the cylinder. 

The crucial property, however, is the fact that the intrinsic equations degenerate on 
the intersections P between the cylinder and null-infinity J !± . Here, the field may 
develop singularities which can propagate onto null-infinity destroying its smoothness. 
Friedrich has conjectured that certain conditions on the initial data have to be satisfied 
in order for the field to be smooth on . 

These analytical properties have been explored numerically in [ 1 1 based on the 
first order system for the spin-2 field. The behavior of the resulting numerical solutions 
near the cylinder I was studied. The performance of the code near the ill-behaved 
location I + at the junction of the cylinder and null-infinity was tested using different 
initial data. The analytical results were reproduced quite successfully. In particular, the 
analytically predicted singular behaviour could be reproduced numerically. The results 
in JTJ strongly indicate that appropriate classes of initial data can be successfully 
evolved without loss of convergence up and including the region I + , but not beyond. 
As expected, since the equations are no longer hyperbolic beyond I + the code develops 
numerical instabilities and crashes. 

With the present article we want to follow up on this work using an equivalent 
system of second order wave equations. Our main motivation for this is the claim by 
several numerical analysts (see e.g., Kreiss and Ortiz |9|) that the numerical solutions 
of second order wave equations have better properties that those for the corresponding 
first order systems. 

The structure of the paper is as follows. In sec. [2] we briefly describe Friedrich's 
representation of space-like infinity of Minkowski space as a cylinder. In sec. [3] we 
derive the system of wave equations and show under what conditions it is equivalent to 
the first order system used in (TJ. In sec.[4]our numerical implementation and results 
are presented. Finally, we conclude with a discussion in sec. [5] The conventions used 
in this work are those of (17). 



2 The finite representation of Minkowski space-time near space-iike infinity 

In this section, we will briefly summarize the finite representation of Minkowski 
space-time close to space-like infinity originally proposed by Friedrich in [6|; more 
detailed discussions of this topic can be found in [ 1 ,6 20). 

The physical metric of Minkowski space-time in Cartesian coordinates (yP ) reads 

8 = Tfyyd/d/, 

where rj^y = diag(l, —1, —1,-1). In this representation the region that we are mainly 
interested in, namely the neighborhood of space-like infinity z , lies far away from the 
origin. In order to get a more detailed description of the points lying at infinity, we 

consider the coordinate inversion Jl6|: y = — which brings our original metric 

f x a x 
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into the form 

Notice that with the aforementioned inversion points close to infinity in our original 
coordinates y 11 are now lying in the vicinity of the origin of the new 'inverted' coordi- 
nates x^ . But this comes at a price: the above metric is singular at the null-cone of the 
origin x^ = 0. Introducing the conformal factor 

Q = -x l x\ 

one can define the conformally equivalent metric 

g' = Q 2 g = ri^dx^dx v , 

which extends smoothly to the null-cone x^x^ = 0. This implies that space-like infinity 
is a regular point in Minkowski space-time. It is well-known that in space-times with 
non-vanishing ADM-mass this is no longer the case and space-like infinity can no 
longer be represented as a point. Instead, as discussed in detail in [6| it is more 
appropriate to represent the region near the origin as a cylinder. This is achieved here 
by a further rescaling of the metric 



with a function fc(r) = rji(r) where r = \J (x 1 ) 2 + (x 2 ) 2 + (x 3 ) 2 is the radial spatial 
distance from the origin and fi(r) is a positive radial function with ji(0) = 1. The 
choice of the function ji represents the remaining conformal gauge freedom. After the 
introduction of a new time-coordinate t by defining x° = K(r)t the metric g expressed 
in polar coordinates takes the following spherically symmetric form 

g = KT 2 (K 2 dt 2 + 2tK fc'dr dr - (1 - 1 2 K l2 )&r 2 - r 2 (df} 2 + sin 2 9 d0 2 )) , (2.1) 

where (0,0) are the usual angular coordinates and ' denotes differentiation with 
respect to the radial coordinate r. Notice that the final rescaled metric d2T) and the 
original Minkowski g are related by the conformal factor 

= K- 1 i2 = -^(l-r 2 M (r) 2 ), 

which is positive on M — {r > 0, |?| < l/jJ,(r)}, corresponding to the Minkowski 
space-time. The conformal factor vanishes on the two regular 3-dimensional hyper- 
surfaces 

J r± = {r>0,t = ±-^r-} and / = {r = 0, \t\ < 1}. 

The hyper-surfaces J !± are null, corresponding to null-infinity, while / is a regular 
hyper-surface with the topology of a cylinder. Interestingly, in the limit r — > future 
and past null-infinity do not meet at the same point as in the conventional picture. 
Instead, they meet the cylinder / in the 2-spheres P = {r = Q,t = ±1}. This obser- 
vation leads naturally to the finite representation of space-like infinity, where 2° has 
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been blown up to the cylinder 7 = {r = 0, — 1 <t < 1}. The function jj. dictates the 
shape of these structures. Here, the two simplest choices will be considered: ji(r) = 1 
and /i(r) = ■A-j.. In the former 'horizontal' representation, null-infinity is given by 
the set J? ± = {r > 0, t = ± 1}, while in the latter 'diagonal' representation by the 
set J± = {r>0,r = ±(l + r)}, see Fig. [T] It is helpful to introduce the double null 




Fig. 1 The neighborhood of the cylinder / for the horizontal and the diagonal conformal representations. 



coordinates 

u = K(r)t — r, v=K(r)t + r, (2.2) 
which puts the metric into the form 

1 1 , 
g = —^dudv ^dar. 

K )1 2 

Now, J !± is characterized by the vanishing of exactly one of the null coordinates, 
the other one being non-zero. Both coordinates vanish on the cylinder I. In Fig. [T] 
we also display the lines of constant u and v. Since these are coordinates adapted 
to the conformal structure of M we can see clearly, that the cylinder is 'invisible' 
from the point of view of the conformal structure. Notice, that the differentiable 
structures defined by the (u,v) and the (f,r) coordinates are completely different near 
the boundary r = 0. 



3 The spin-2 wave system 

Our objectives in the present section are to derive a second order system for the 
spin-2 equation described in JTJ and to establish a correspondence between the two 
formulations. 
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3.1 The spin-2 wave equation 

Our starting point will be the spin-2 zero-rest-mass equation (16) on an arbitrary 
space-time 

V E A«l>EBCD = 0. (3.1) 

Acting on ((XT) with another spinorial covariant derivative V^a' yields 

- V^,V £A > £BCD = 1 -n<j> ABCD + 3y A{B EF <j> CD]EF +6A<j> ABCD . (3.2) 

In the present case, the Weyl spinor vanishes, since the metric (2.1) is conformally flat. 
The Ricci tensor of this metric turns out to be trace-free, so that also A = 0. Thus, we 
obtain the wave equation 

a(j) ABCD =0. (3.3) 



3.2 Relating the first with the second order system 

Now, we will relate the first ( |3.1) with the second ( |3.3[ ) order system. Establishing 
this correspondence is necessary because we need to know under what conditions we 
obtain the same solutions. We will then also be able to compare the numerical behavior 
of the two systems and study their similarities and differences. In the following, we will 
prove that in a conformally flat space-time with vanishing Ricci scalar a solution <p AB cD 
of the spin-2 wave equation ( |3.3| >, which satisfies initially the spin-2 zero-rest-mass 
equation pTTT l, is also a solution of the latter. 
First, we define the spinor 

Za'BCD = V F A> §FBCD ■ (3 .4) 

Then the equation E A i BCD = is nothing else than the spin-2 equation p.l) . Rewriting 
( |3.2) , we obtain an equation for the derivative of E A < BCD 

Vaa'^ A 'bcd = \a$ ABCD + 3>¥ A{B EF 4> CD)EF +6A$ ABCD . (3.5) 
Now, assume that (j) AB co satisfies (3.1) , then E A / BCD = and 

hj(j> AB CD + ^A(B EF ^CD)EF + 6A<j> ABC D = 0. 

Contracting over indices A and B yields 

V(C AEF $D)AEF=0, 

a purely algebraic condition between the conformal curvature and the spinor field. This 
is a Buchdahl condition fTT) which implies the well-known fact that the zero-rest mass 
equation (with spin 2 in this case) is inconsistent on space-times with non-vanishing 
Weyl tensor. So we take *¥ AB cd = from now on and we also assume that A = 0. Then 
the spinor field satisfies the second order wave equation 



n<j) ABC D = o. 



(3.6) 
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Conversely, let §abcd be a solution of ( |3.3| l. Then ( |3.5| ) shows that the spinor field 
La' bcd — satisfies the equation 

V^Z A 'bcd = 0. (3.7) 

In order to study the properties of this system of PDE we look at its symbol. This 
is, for every real co-vector paa', defined as a map P : y A ® -5^(bcd) ~^ ^a ® ^{bcd) 
between the appropriate spin spaces defined by 

a A 'PBCD i-> Paa>u A 'Pbcd- 
Therefore, it is enough to discuss only the map 

a A ' ^ pAA/a A ' . (3.8) 

Choosing a spin basis (o A ' , l A ' ) in y A> and the corresponding basis (— l A , oa) m -^A 
and expanding p^A' as 

A4A' = Po (o^oa' + UU') - P\ {oa Ia' + oa U' )-&p2{o A l A i-o A l A i)-p?, {o A o A i - l A l A , ) 
with real po,. .., pi we find the matrix representation of the map ( |3.8| > 

( PO + Pi Pi -&P2\ 
\PI+Rp2 po-pi J 

which is clearly hermitian for every choice of co-vector paa' an d it is positive definite 
for the co- vector with p\ — pi — pi — 0, po — 1. The matrix representation of the 
map P then simply consists of four copies of this matrix along the diagonal. Thus, the 
system is symmetric hyperbolic and its Cauchy problem is well-posed. 

Thus, if Ea'bcd = on an initial hyper-surface and, for initial boundary value 
problems, if boundary data are given such that the ingoing Ea'bcd waves vanish, then 
the field Ea'bcd vanishes everywhere, i.e., the field Qabcd satisfies pTTJ. This conclu- 
sion holds on the region where the systems are hyperbolic. In our application below 
we will consider situations where hyperbolicity is lost on a region with codimension 2, 
namely the spheres / . We find empirically that the conclusion still holds near these 
regions. 

We will use this result as follows. In the numerical implementation of ( |3.3[ ) the 
initial data and the boundary conditions are determined from ([XT}. Then, the computed 
solution will be compared to a previously computed |[TJ solution of < |3 - 1 1 > - 



3.3 Coordinate representation of the second order wave equation 



In order to obtain a coordinate representation of ( |3.3| ), we introduce a null-tetrad that 
is adapted to the spherical symmetry of the background metric ( |2.1| >. We will choose 
the same null tetrad (Z v ,n v ,m v ,m v ) as in |1J, namely 

V = -^(l-ffc',K-,0,0), n v = -^(l+fK-',-K-,0,0), 
v 2 v 2 



M o,o,i,-^ 1 - M 



^2 V ' sine/' V2 V 'sin0 



0,0,1, 
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The non-vanishing spin-coefficients for this tetrad are 

p = -p' = ~r M ', e = r =~K'. (3.9) 

Now we express the wave operator □ in terms of the weighted differential operators 
of the GHP (Geroch-Held-Penrose) formalism fT7) using the commutation relations 

(p'p - pp')0ABCD = 0, (3p - p3)0ABC£> = ~P 

(W -VWABCD = -p'&faBCD, (S'3 - W^ABCD = o, 

which contain only the non-vanishing spin-coefficients. Expanding the spin-2 zero- 
rest-mass field in the familiar way 

<!>ABCD = lAlfilcl£)0O — 4l( A l B lcO D )(^l + 

+ 61( A 1bo c o d) 2 -4i (a o b o c o d) 03 +o a o b o c o d $ a , 
we obtain the following system of five equations 
W^k-p'Mk-p^'<l>k+k(5~k)p 2 <j) k ^dd'<j) k + (4-k)pd<j) k+l ~kpd'<j) k ^ (3.10) 
with k — 0, 1,2,3,4. 

Expressing the GHP operators p and p' in terms of the directional derivatives along 
the tetrad vectors I and n (see [ 17 1) yields 

pT] = - 1 ={{l-tK , )d t + Kd r -2\/2we)T}, 
V2 

p'rj = -^-={{\+tK')d t -Kd r -2V2wy)7}, 
v2 

where r\ is a {/^gj-scalar quantity with boost-weight w — and spin-weight 
s = ^2^, see Jl7j. Due to the spherical symmetry of the metric and the adapted 
null-tetrad it is natural to refer the 3 and 3' operators to the unit-sphere, so that we 
replace 

g^Ag, g'^Ag' 
V2 V2 

in ( |3.10| ) as in Q. Then the system ( |3.10| ) takes the form 

{l-t 2 K a )d tt <$) k -K 1 d rr <$> k + 2tKK'dtr<t>k 

+ 2[{2-k)x' -t^ + r/lV -^KK")]d t <j> k + 2rKn'd r ^ k 

+ [(2 - k){KK" + (l-k) K 12 ) + Jfe(5 - k)r 2 pL 12 ] fa 

= H 2 d& <j) k + (4 - k) r /J. n' B(f) k+1 - kr n n ! & <j) k _x. 

Finally, we use the spherical symmetry of the metric ( |2.1| > to expand the components § k 
of the spin-2 field as a sum of spin-weighted spherical harmonics s Y[ m in the following 
way 

#(f,r,0,0)=^^ m (f,r) 2 _,y /m (0,0), 

llll 
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where s = 2 — k is the spin-weight of fa an d me integers s,/,m satisfy the inequali- 
ties |s| < I and \m\ < I. Since the operators 3,3' act on the spin-weighted spherical 
harmonics s Y[ m as 

%Y lm ) = -y/lQ + l)- s ( s+ l) s+1 Yi m , 
&( s Y lm ) = y/l(l + l)-s(s-l) s -iYi m , 
SSVim) = -(1(1 + 1) S Ylm, 

the system ( |3. 1 1 ] > decouples into separate systems for each mode of the fixed pair 
(l,m), i.e. 

(1 - t 2 K l2 )d t ,fa - K 2 d rr fa + 2t KK'dtrfa 

+ 2[(2-k) K ' -t(K a + r/iV - l -KK")}d t fa + 2rKn'd r fa 

2 (3. 12) 

+ [(2 - k)(KK" + (1 - Jfc) K"' 2 ) + ifc(5 - k)r 2 pL l2 ]fa 

= -n 2 c 2 k fa- r^n '((4-k)c k fa +1 +kc k ^fa^i), 

where c/, = ^Jl(l + 1) — (2 — k)( \ — k) whenever the square root is real and otherwise 
q = 0. For the sake of simplicity the notation = fa was introduced. It is worthwhile 
to mention that this form of the system holds for / > 2, since for / = 1 the components 
00 and fa vanish and for / = only fa survives. Since the above system comes from 
equations involving the wave operator it is hyperbolic. 



3.4 Characteristic curves 

Because of its hyperbolic nature the system ( |3. 12[ > has two families of real characteris- 
tic curves. As in the case of the spin-2 equation, their behavior will be very useful in 
the subsequent numerical studies. Following |8], the slope of the characteristic curves 
for the system ( |3.12[ ) is given by 

dt l-tx!(r) 
dr K(r) 

The characteristic curves of the spin-2 wave equation are identical to those of the 
spin-2 equation. Fig.[2]shows the 'outgoing' characteristic curves of ( |3. 12) > with slope 
dt /dr = (1 — tK 1 ) I K near a neighborhood of 7 + ; near I these curves tend to become 
parallel to the cylinder I. On the other hand, the 'ingoing' characteristic curves with 
slope dt /dr = — (1 +t k')/k are becoming parallel to I near 7 + ; their behavior in a 
neighborhood of I can be visualized by a reflection of Fig. [2] along the r-axis. 

There is a slight difference with the first order case, though. Now, the behavior of 
the characteristic curves is universal, in the sense that all the components of the spin-2 
field propagate along the same characteristics, while in the first order system different 
components have in general different characteristics (depending on the gauge). 

Furthermore, as in the first order system the evolution equation reduce to an 
interior system on the cylinder in the sense that there remain no r-derivative terms in 
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Fig. 2 Characteristic curves of all the fields ^ in a neighborhood of / + . The red line denotes the cylinder 
and null-infinity. The situation close to I~ is obtained by a reflection in the r-axis. (a) Characteristic curves 
of slope = ^y- in the 'horizontal' case fl = 1. (b) Characteristic curves of slope jjc = 1- ' in the 
'diagonal' case jl = . 



the equations when evaluated at r = 0. So, as before, the cylinder is a total characteristic 
for the system. And as before, the equations loose their hyperbolicity at / , since the 
coefficient in front of the second time derivative vanishes there. 



4 Numerical analysis and results 



In this section, we will numerically solve the initial boundary value problem for the 
spin-2 wave equation ( |3. 12[ > using the findings of sec. 3.2 Specifically, we will pre- 
scribe initial data that satisfy the constraints of the first order system and subsequently 
evolve them with ( |3. 12[ >. In addition, the boundary conditions will be prescribed using 
the first order system. 



4. 1 Numerical implementation 

We discretize the 1 + 1 system ( |3.12| > (one time and one spatial dimension) using the 
method of lines. Accordingly, the system ( |3.12| > is reduced to a system of ordinary 
differential equations by discretizing the spatial coordinate r. In this work, we use 
finite difference techniques to achieve that. 

As computational domain we adopt the interval D = [0, 1]. A finite representation 
of D is obtained through the introduction of an equidistant grid r, = ih with i = 
0, . . . ,N onD, where r# = 1, and h is the grid spacing. The spatial derivative operators 
appearing in ( |3.12| i are approximated by appropriate summation by parts (SBP) finite 
difference operators [3.4, 15 19 J. The SBP operators arise from discrete versions of the 
continuous energy estimates that were originally introduced in [10, 11]- While these 
energy estimates guarantee the well-posedness in the continuous case, on the discrete 
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level their discrete counterparts guarantee the numerical stability of our discrete 
schemes. This is the main reason we chose to use these specific finite difference 
operators. 

Another very appealing property of the SBP operators is that although their 
accuracy near the boundaries is, depending on the details of their construction, one 
or two orders smaller than the one in the interior of the grid, the overall accuracy is 
of the same order with the accuracy in the interior. Here, we approximate the second 
derivatives with a minimal width, full norm SBP operator given in [15], which is 
fourth order accurate in the interior and second order at the boundary. Unfortunately, it 
was not possible to stabilize the above operator with the corresponding SBP operator 
of the same norm matrix approximating the first derivatives given in |T3] . Instead, we 
combined it with the SBP operator constructed in fl9) , which was very successfully 
used already in the numerical implementation of the first-order formulation of the 
spin-2 equation, see | Thus, the first derivatives were approximated with a minimum 
bandwidth, restricted full norm SBP operator of fourth and third order accuracy in the 
interior and at the boundaries, respectively, taken from 1 19 1. 

With this choice, we violate the assumption made in p3] , of approximating both 
first and second derivatives with SBP operators of the same norm. Theoretically 
though, under certain conditions on the norm matrices of the two operators, energy 
estimates like the one at the end of sec. 14.21 can still be obtained. So it seems that 
requiring the same norm for first and second order operators is sufficient but not 
necessary for stability. To thwart any possible objections — emanating from the above 
complication — concerning the performance of our numerical code we also carried out 
numerical simulations using (first and second order) SBP operators of the same norm — 
but one order less accurate than the choice above. Specifically, the second derivatives 
were approximated by a minimal width, diagonal norm, fourth order accurate in the 
interior and second order at the boundary SBP operator and the first derivatives by a 
SBP operator of the same accuracy and norm, both given in p3| . Being constructed 
for a diagonal norm this combination of operators is expected (see [ 15 1) to have a third 
order overall accuracy. Our numerical findings, see |5}, confirmed this expectation. 

A point that also needs special attention is the imposition of the boundary con- 
ditions. A wrong imposition of the boundary condition would destroy the designed 
accuracy of the SBP operators and lead to instabilities [2|. We use a very simple, 
but highly efficient, penalty method — the so-called simultaneous approximation term 
(SAT) method introduced in |2} — that preserves the designed accuracy of the SBP 
operators and guarantees the numerical stability of our schemes. The combination 
of the SBP operators with the SAT method has been successfully applied in several 
different circumstances, see e.g., pf|4][7l [T2l[T5l[T8) . 

Finally, one has to decide how to solve the semi -discrete system of ODEs obtained 
from p,12[ ) after the spatial discretization. Here, we implement the system as a first- 
order system of ODEs by introducing the first order time derivatives of the grid 
functions as additional variables. This system is then solved using the standard 4th 
order Runge-Kutta scheme. 

We choose to work with the first-order system for a number of reasons. Firstly, its 
numerical implementation is quite simple and very well studied. The imposition of the 
boundary conditions with the SAT method is also simpler in that case. Furthermore, 
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a second-order system of ODEs imposes additional restrictions on the construction 
of SBP operators, which limits the operators available to those with diagonal norm, 
see JT5) for the details. 



4.2 Treatment of the boundaries 



As it was already mentioned in sec. 3.4 the cylinder is a total characteristic of the 
system ( |3.12| i. Thus, we do not have to prescribe boundary conditions for the points 
of our computational domain that lie on the cylinder, i.e. at r = 0. In contrast, at the 
boundary r = 1 there is one ingoing characteristic for each component of the spin-2 
field. Thus, we have to provide appropriate boundary conditions at r = 1 for each 
component This is an important difference to the first-order spin-2 equation ( |3.1[ ), 
where only some of the components propagate inwards and therefore need boundary 
conditions. We use different boundary conditions which are known to lead to well- 
posed problems for the wave equation: Neumann or, more generally, Robin conditions 
of the form 

<l> k (t,l) + d r <l> k (t,l)=g k (t), (4.1) 

where k = 0, . . . ,4, as well as Sommerfeld type boundary conditions. Since we need to 
use informations from the 1st order formulation of the spin-2 equation (as discussed 
in sec. |3.2[ l, which provides relations for the values of the first derivatives of the 
components but not for the values of the functions themselves, we cannot impose 
Dirichlet conditions. 

To motivate the above choice of boundary conditions and to exemplify the use of 
the SAT method for systems with second derivatives, we will consider the simplest 
form of a 1-D wave equation 

u„ =u xx , < x < 1 , t > (4.2) 

with homogeneous Robin boundary conditions 

au(t,0) + u x {t,0) = 0, J3k(?,1) + 1) = 0, 

where a,/3 are arbitrary constants. The energy method for the above boundary value 
problem gives 



±(j\u 2 + u 2 x )dx^ =j t (-pu 2 (t,l) + au 2 (t,0)), 



which leads to an energy estimate when j3 > and a < 0. 

The semi-discrete approximation of eq. ( |4.2| i is i)„ = D2V, where v — (Uj)i=o:iV is 
the grid function approximating the solution u and the (N +1) x (N +1) matrix D2 
is the SBP operator approximating the second derivative operator d 2 /dx 2 . Following 
(T3"][T5) , the implementation of the boundary conditions with the SAT method leads to 
the introduction of a couple of boundary terms in the following fashion 

v tt = D 2 v + To H - 1 E ( a I + S ) v + t n H - 1 E N ( j3 / + S ) v , 
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where Eq = diag(l,0, . . . ,0), En — diag(0, ... ,0, 1), / = diag(l, . . . , 1) are matrices 
of size (N + 1) x (N + 1 ), S is an approximation of the first derivative operator at the 
boundary, H is the norm matrix of D 2 , and To, % are the so-called penalty parameters. 
The energy method for the above semi-discrete wave equation gives 1 13 14) 

j (v?Hv, + v T Av) = j t {z avl ) + T N j5v^)+2(t Q -l)vo{Sv)o+2{z N + l)v N {Sv) N , 

where A = A T is the matrix in the definition of D 2 , i.e. D 2 — H~ l (—A +BS) with 
B = diag(— 1,0,. ..,0,1). An energy estimate exists for To = 1 , In — — 1 an d (as above) 
P >0, a<0: 

j t {vjHv, + v T Av) = j t {avl-$v 2 N ), 

which is the semi-discrete analog of the above continuous energy estimate. This 
procedure for imposing boundary conditions holds true [ 14 1 also for Neumann and 
Sommerfeld conditions, but not for Dirichlet conditions. In that case, as was shown 
in |14| , the method of imposing boundary conditions with the SAT method must be 
modified in order to guarantee the stability of the numerical scheme. As discussed 
above, we do not have the need for imposing Dirichlet conditions. Therefore, by 
choosing ( |4.1| i, Neumann or Sommerfeld boundary conditions, we do not have to go 
into these additional complications. 



4.3 The exact solution 

As described in sec. [3] the system ( |3.12[ ) is just a second order reformulation of the 
spin-2 equation. Therefore, it must be satisfied by the exact solution presented in (T]|5): 

. , , r 2 (l+r-f) 4 . , , 2r 2 (l+r-f) 3 (l+r + f) 
*>(',') = (1+r)7 , M,r) ^ , 

V6r 2 fl +r-t) 2 (l +r + t) 2 
(1+r) 7 

2r 2 (l+r-f)(l+r + f) 3 r 2 (l+r + t) 4 



<h(t,r) = , fa(t,r) 



(1+r) 7 ' r " v ' ; (1+r) 7 • 

Recall that this solution refers to the 1 = 2 mode in the 'diagonal' representation of 
null-infinity. According to the discussion in sec. 3.2 the initial data must satisfy the 
constraints of the first order system and be subsequently evolved with the second order 
system ( |3.12| >. Therefore, evaluating ( |4.3[ ) at t = 0, the initial data for each one of the 
fields are 

^(r) = 4 (r) = ^^3, fr( r )=k( r ) = ^_ fc ( r) = _|L_. (4 . 4 ) 



In a similar way, i.e. differentiating ( |4.3[ ) with respect to t and evaluating at t = 0, the 
first temporal derivative of the components of the spin-2 field read 

4r 2 

(iT7) 4 



Mr)=fa{r) = -k{r) = -h{r) = - T — ?) <j> 2 (r)=0. (4.5) 
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In summary, ( |4.4| i and ( |4.5| l will be our initial conditions in this section. In addition, 
at r = 1, we have to prescribe Robin boundary conditions of the form ( |4. By 
differentiating ( |4.3[ ) with respect to the spatial coordinate r and evaluating at r = 1 we 
obtain 

8o{t) = " ?k {t ~ 2)3 (6 + 1] ' gl (f } = m (48 _ 32f + f4)> 

. 2 W = -4/|(^-4)(12 + ^), 

53(0 - (48 + 32* + » 4 ), 84(t) - - ^(r-6)(2 + f) 3 , 

where the subscript denotes the component of the spin-2 field to which each boundary 
condition corresponds. The above boundary conditions are imposed with the SAT 
method. 





00 


04 


Grid 




1 st order 


2nd order 


1st order 


2nd order 


50 




-24.8993 


-27.1097 


-11.3229 


-12.2922 


100 




-29.0579 


-31.3535 


-13.3223 


-15.1472 


200 




-33.3784 


-35.5212 


-17.3260 


-18.1698 


400 




-37.7813 


-39.6182 


-20.4876 


-21.2892 



Table 1 The logarithm of the normalized I 2 norm of the absolute error £, log 2 (| \E\ 2), between the exact 
solution and the solutions computed from the lst-order system (3.1) and the 2nd-order system {33) at time 
t = 1. Note, the 4th order convergence in 0o and the 3rd-order convergence in 04. 



Evolving the initial data ( |4.4| > and ( |4.5| l with the second order system ( |3.12| > we 
can reach t = 1 without loss of the expected 4th order convergence. In addition, the 
constraints are preserved during the evolution. The above results strongly indicate that 
our code reproduces successfully the exact solution ( |4.3| l. 

But the purpose of this work is to compare the present numerical approach with 
the one developed in [1], thus a first step towards this goal would be to check which 
approach reproduces better the exact solution ( |4.3| l. Tab.[T]serves this purpose. Specifi- 
cally, it depicts the logarithm of the normalized I 2 norm of the absolute error E, i.e. 
log 2 (| \E\ I2), of the components 0o, 04 for the two numerical approaches at time t = 1. 
Notice that a slightly better accuracy is achieved in the second order approach. 

We have also done the comparisons in the 'horizontal' case with similar results. 
Here, we cannot reach I + due to the degeneracy of the equation at t = 1 which affects 
the numerical algorithm because the propagation speeds grow unboundedly. However, 
we can come arbitrarily closely to using an adaptive time-step. We will not go 
into any further details here because the phenomena are the same as discussed in |TJ. 



4.4 General initial data 



According to the results of the preceding section, there is strong evidence that our 
numerical code converges with the expected order and reproduces successfully the 
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exact solution ( |4.3| l of the system ( |3.12| i. With the confidence that these results provide 
us, we can proceed further in the numerical study of the spin-2 wave system and 
seek for numerical solutions, which do not correspond to exact solutions of the 
system ( |3.12| >. In addition, we will confirm numerically the analytic result of sec. 3.2 



namely that the solutions of the two different formulations for the spin-2 equation 
coincide if the initial and boundary data are chosen appropriately. 

Again, the initial data satisfy the constraints of the 1st order system and are evolved 
with the second order system ( |3.12[ ). We use the constraints in the form derived in 
sec. 3.3 of fl] to determine the initial values of all the components in terms of one 
free function. We define the auxiliary variables y/,(r) := 0,-(r)/jU(r) 3 for which the 
constraints imply the following relationships 



Vb(r) = Mr) 



- a$y 2 W + 2 r (r) + 2 r 2 \j/ 2 ' 



aoa 2 



•v4( ? 

do 



(4.6) 



with a = y/lQ + l) and a 2 = ^l(l+l)-2. 

This choice leaves the component 02 completely at our disposal and allows us to 
compute the remaining components explicitly. Note, that this is not the most general 
form of the initial data. We have restricted ourselves to the case where 03 (r) = 0i (r). 
Otherwise there would exist another free function. 

We take the initial values of 02 in the form of a bump function 



ur)= \T bK b,) ' , (4-7) 




centered &tr = b/2, then all the other components are bump functions as well. In this 
section, we will choose b = 1. 

We also have to specify initially the first time derivatives of the components. The 
evolution equations of the first order system, see Appendix |A| will be used for this 
purpose. Evaluated at t = they read 

00 (r) = K0o( r ) - ( 3k> -M) <h( r ) ~ °<2M0i i r ), 

02(r) = -ab|U0i(r)-iao/i<Mr), ( 4 -8) 
fa{r) = -CColl<l>2(r)-^a2ll<l>4(r)+ll<l>3(r), 



where the values of the fields on the r.h.s can be evaluated from ( |4.6| l, ( |4.7| i. 

Finally, we explain how we specify the boundary conditions. As discussed in 
sec. 3.2 we need to use the available information from the first-order system ( |3.1[ ). In 



this system there is only one component, namely 0o, which propagates inward from 
the boundary at r = 1 . Thus, there is only one free function to be specified on that 
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boundary, which characterizes the solution inside given initial data. This must be the 
case also for the second order system after imposing the boundary conditions. 

The second order wave equations require for each component a boundary condition 
at r = 1 (recall that r = is a total characteristic so we cannot prescribe any conditions 
there). We impose these conditions in the form of a Robin condition ( |4. l| i 

fc(f,l) + #(t,l)=ft(t)» i = 0,1,2,3,4. 

The boundary functions gj(t) are computed from the first-order system in terms of the 
fields and their time derivatives. There is no unique way to do this because there are 
both constraint and evolution equations involving 0i, 02 an d 03 which could be use 
for this purpose. We choose to use the constraint equations ( A.6[|A.8 1 for 0i, 02 and 
03 and the evolution equations dAT) and (|A.5[) for 0o and 04. 



Since 0o is freely specifiable on the boundary we choose it simply as zero, so 
that 0o and its time derivative vanish on the boundary. Wherever 0o appears in the 
equations used, we simply drop it. The functions gj(t) are computed by evaluating the 
left hand side of ( |4.1[ ) with the spatial derivative substituted from the corresponding 
first-order equation and putting 0o and 0o equal to zero. 

The resulting equation for each component is imposed numerically using the SAT 
method. Note, that we have essentially only rewritten the first-order equation for the 
components in a superficial form of a Robin condition. We could just as well have 
chosen a Neumann condition or even a Sommerfeld condition. In all cases the effective 
boundary condition, i.e., the penalty term to be added to the discretized wave equation 
would have been the same. In fact, we do not see any difference if we implement the 
same boundary conditions as Neumann or even as Sommerfeld condition. 

We specify initial data for the simplest non-trivial case / = 2 in the 'diagonal' 
representation. For the evolution of the data ( |4.7| i, ( |4.8[ ) described above, the numerical 
solutions for the 'ingoing' and 'outgoing' components 0o and 04, respectively, are 
presented in Fig. [3] These components satisfy first-order advection equations which 
are purely ingoing resp. outgoing. These properties are clearly visible in the plots. 
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Fig. 4 The convergence plots of <j>4 for the evolution of the initial data (4.7) , (4.8) in the 'diagonal' 
representation at time t = 1 in (a) the first order approach [1] and (b) the second order approach developed 
in the present work. 



In Fig. [4] we compare the convergence plots of the 04 component at t = 1 obtained 
from the evolution of the same initial data ( |4.7| i and ( |4.8| l in the two numerical 
approaches developed in the present work and [ 1 1. More precisely, in Fig. 
initial data are evolved using the first-order system (|A. 1]>-([A75]>, while in Fig. 



4(a) 



4(b) 



the 
the 

second order system ( |3. 12) > is used. In each case, the difference to a high resolution run 
is computed, which here is 800 grid points. In both cases we find 4th order convergence. 
Notice, though, that in the second order case the plots look much cleaner. The high 
frequency features, appearing on the convergence plot of the first order system, are 
not present in the plot for the second order system. This result confirms the prediction 
of (|9) that the spurious waves will disappear in the second-order formulation. 





<k 


<j> 4 


Grid 


log 2 (||£|| 2 ) Rate 


log 2 ( \\E 2) Rate 


50 


-2.5929 


-1.2116 


100 


-6.4693 3.8764 


-5.1393 3.9277 


200 


-10.3956 3.9263 


-9.1361 3.9968 


400 


-13.8513 3.4557 


-13.3382 4.2021 



Table 2 Convergence of the solutions to the wave equations towards a solution at high resolution (800 
grid points) of the first order equation. The table shows the (logarithms of) the normalized I 2 norm of the 
absolute differences and the corresponding convergence rates at time f = 1 . 



Now, we reproduce numerically the analytic result of sec. 3.2 we show that if 



the initial data satisfy the constraints ( |4.6[ > of the first order system and are evolved 
with the second order system ( |3.12[ ), then the resulting numerical solutions converge 
(with increasing resolution) to the one obtained by evolving the same initial data with 
the evolution equations ( |A. l| i-( [A75t > of the first order system. To make a connection 
with the results of Fig.|4j we consider again initial data centered around r — 0.5 of 
the form \A.1\ , ( |4.8| l and impose boundary conditions as described above. We evolve 
these data with both evolution systems and compute the absolute error between the 
numerical solutions of the second order approach and the numerical solution of the 



The second order spin-2 system in flat space near space-like and null-infinity 



17 



highest resolution (800 grid points) in the first order approach. The normalized I 2 
norm of the computed absolute error and the corresponding convergence rates at time 
t = 1 for the components 0o> 04 are listed in Tab. [2] We find that the solutions agree 
within numerical accuracy. This confirms the statement made in sec. 3.2 that the two 
systems have identical solutions given the same initial data. 



5 Discussion 

In this work we have reformulated the first-order system for the spin-2 field equations 
for linear gravitational perturbations on Minkowski space as a system of coupled 
second-order wave equations. The system was implemented into a numerical code and 
studied under certain simplifying assumptions. 



Our first analytical result, see sec. 3.2 is that the first-order system and the second- 
order system are equivalent if and only if the underlying space-time is conformally 
flat and has a vanishing scalar curvature. In the case under study the space-time is 
conformal to Minkowski space-time, i.e., conformally flat and the conformal factor 
has been chosen in a such a way that the scalar curvature vanishes. Thus, we can 
evolve the spin-2 field using either the first or the second order system. Given initial 
and boundary data computed from the first-order system the second-order system 
provides the same solutions. 

Our main goal was the comparison of the numerical properties of the two formula- 
tions in the conformal setting given in sec. [2] From the numerical results presented in 
sections 4.3 and 4.4 we can conclude that we get the same qualitative behaviour of 
the numerical solutions. In particular, we have shown explicitly that solutions of the 
second-order system converge to a solution of the first-order system with the same 
initial and boundary data. The solutions can be evolved stably up to and including I + 
in the diagonal case but not beyond. Any attempt to evolve beyond I + with the same 
algorithm end in numerical instability. In an upcoming paper we will discuss some 
attempts to circumvent these problems. 

Concerning the comparison of the two approaches we confirm the claims in Q: 
with the same number of grid points we obtain better accuracy roughly by a factor 4 
with the second-order system and we do not see any spurious high-frequency waves in 
the second-order system. 

In our code we have chosen a combination of first and second order discrete opera- 
tors which are SBP operators with respect to different inner products. Our numerical 
results indicate though that the requirement of using consistent SBP operators that 
are based on the same norm does not seem to be very restrictive. The fact that we 
could obtain results of similar accuracy and convergence with the ones presented in 
secs. |4.3||4.4| with yet another restricted full norm operator taken from [4] strengthens 
this statement even more. It seems that the second order full norm operator we are 
using is performing quite well when combined with first order SBP operators based on 
a restricted full norm. We have also carried out numerical simulations using (first and 
second order) SBP operators of the same norm. Our results agree with their designed 
accuracy as described in 1 15 1. 
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A The first order equations 

As a reference, we list here the spin-2 zero-rest-mass equatio n <|3 . 1 ^ in the formulation used in fTJ. The 
geometry and the notation is exactly the same as described in sec.|3.3| The system consists of eight equations 



which can be split into the five time evolution equations 

(l+fK J )d,fo-lcd r 0o = -(3K'-M)0O-Mra20l! (A.l) 

d,tj>i = + ^a 2 (j>o- ^/x«o02, (A.2) 

d,tpi = ^Aiao0i - ^/iOb03. (A3) 

3,f 1 = ^03 + ^a o <|)2-i/xa2^4, (A.4) 

(1 -tK')d,<p4 + Kd r <j>4 = (3ie'-/i)04 + tia 2 <fe (A.5) 

and the three constraint equations 

-2Kd r <pi +6r/i'0i -ItK'lifa +Oo/l(l -ffc')02 + O!2/x(l +/k')0o = 0, (A.6) 

-2Kd r (j> 2 +6rii'<l>2 + a ii(i-tK')(j> 3 + a o ii(\+tK')(j> [ = 0, (A.7) 

-2Kd r fo +6r/i'0 3 +2/K / /i</> 3 + Oo/l(l + ffc')02 + a2/i(l -tn!)fa = 0. (A.8) 
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